#########################################################################
# Author:	Olga Chyzh and Frederick J. Boehmke		 	#
# Date:		May 14, 2015				   	  	#
# File:		table1-ergm.do					  	#
# Purpose:	Replicate Cranmer and Desmarias' (2011) replication of	# 
#		Long (2008). Note that the results with the imputed	#
#		data are outputted to be read in and combined later.	#
# Input File:	table1-ergm.dta						#
# Output File:	table1-ergm-base.txt,					#
#		table1-ergm-ivnetworkX.txt,				#
# 		table1-ergm-ivnetworkX.txt,				#
# 		table1-ergm-ivscoreX.txt				#
# Previous:	This was adapted from Cranmer and Desmarais' 		#
#		replication for their article.				#
#########################################################################

# Change to the working directory

setwd("D:/International/confnetworks/journals/psrm/replication files")

# Load the needed libraries.

library(foreign)
library(sna)

# Declare the variables and data sets.

maoz.vars <- c("year","statea","stateb","dichmid20","Lminreg302","Lcaprat301","Ldistance","Latopsecorr","Ligosecorr","Ltradesecorr","Lintegcorrse","midspline1","midspline2","midspline3","t1","t2","t3","Ltrade_yhat_seq1","Ltrade_yhat_seq2","Ltrade_yhat_seq3","Ltrade_yhat_seq4","Ltrade_yhat_seq5","Ltrade_seq_hat1","Ltrade_seq_hat2","Ltrade_seq_hat3","Ltrade_seq_hat4","Ltrade_seq_hat5","peaceyrs")

  maoz <- read.dta("table1.dta")
  maoz.dat <- NULL
  for(i in maoz.vars){
	maoz.dat <- cbind(maoz.dat,maoz[,which(names(maoz)==i)])
	}

  maoz.dat <- data.frame(maoz.dat)
  names(maoz.dat) <- maoz.vars

  maoz.dat <- na.omit(maoz.dat)

# Create triadic changestats

tchange <- numeric(nrow(maoz.dat))
tschange <- numeric(nrow(maoz.dat))

for(i in 1:nrow(maoz.dat)){
	dat1 <- subset(maoz.dat,maoz.dat$year==maoz.dat$year[i])
	dat1mid <- subset(dat1,dat1$dichmid20==1)
	datc1 <- subset(dat1mid,as.numeric(dat1mid$statea==maoz.dat$statea[i])+as.numeric(dat1mid$stateb==maoz.dat$statea[i]) > 0)
	datc2 <- subset(dat1mid,as.numeric(dat1mid$statea==maoz.dat$stateb[i])+as.numeric(dat1mid$stateb==maoz.dat$stateb[i]) > 0)
	us1 <- unique(c(datc1$statea,datc1$stateb))
	if(is.element(maoz.dat$statea[i],us1)) us1 <- us1[-which(us1==maoz.dat$statea[i])]
	if(is.element(maoz.dat$stateb[i],us1)) us1 <- us1[-which(us1==maoz.dat$stateb[i])]
	us2 <- unique(c(datc2$statea,datc2$stateb))
	if(is.element(maoz.dat$stateb[i],us2)) us2 <- us2[-which(us2==maoz.dat$stateb[i])]
	if(is.element(maoz.dat$statea[i],us2)) us2 <- us2[-which(us2==maoz.dat$statea[i])]
	if(min(c(length(us1),length(us2)))>0){
		tchange[i] <- length(intersect(us1,us2))
		tschange[i] <- length(us1)+length(us2)
	}
}

whichab <- function(x,y){

  return(which(y==x))

  }

# Bootstrap year-by-year.

maoz.dat <- cbind(maoz.dat,tchange,tschange)

  m <- 1000
  set.seed(10)
  uyr <- unique(maoz.dat$year)

# Run the base model without instruments to compare to previous results and our modifications.

coef.triad0 <- NULL
coef.triad1 <- NULL
coef.triad2 <- NULL
coef.triad3 <- NULL
coef.triad4 <- NULL
coef.triad5 <- NULL

for(i in 1:1000){
	yrsi <- sample(uyr,length(uyr),rep=T)
	yrsil <- as.list(yrsi)
	inds <- unlist(lapply(yrsil,whichab,y=maoz.dat$year))
	datai <- maoz.dat[inds,]
	mod.triad0 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltradesecorr+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)
	coef.triad0 <- rbind(coef.triad0,rbind(mod.triad0$coefficients))
	}

# Output the results to access and combine externally.
	
  dput(list(coef=coef.triad0),"table1-ergm-base.txt")

# Run the models with the instrumented network.

coef.triad1 <- NULL
coef.triad2 <- NULL
coef.triad3 <- NULL
coef.triad4 <- NULL
coef.triad5 <- NULL

m <- 1000
set.seed(10)
uyr <- unique(maoz.dat$year)
for(i in 1:1000){
  
	yrsi <- sample(uyr,length(uyr),rep=T)
	yrsil <- as.list(yrsi)
	inds <- unlist(lapply(yrsil,whichab,y=maoz.dat$year))
	datai <- maoz.dat[inds,]

	mod.triad1 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltrade_yhat_seq1+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)
	mod.triad2 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltrade_yhat_seq2+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)
	mod.triad3 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltrade_yhat_seq3+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)
	mod.triad4 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltrade_yhat_seq4+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)
	mod.triad5 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltrade_yhat_seq5+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)

	coef.triad1 <- rbind(coef.triad1,rbind(mod.triad1$coefficients))
	coef.triad2 <- rbind(coef.triad2,rbind(mod.triad2$coefficients))
	coef.triad3 <- rbind(coef.triad3,rbind(mod.triad3$coefficients))
	coef.triad4 <- rbind(coef.triad4,rbind(mod.triad4$coefficients))
	coef.triad5 <- rbind(coef.triad5,rbind(mod.triad5$coefficients))

	}

# Output the results to access and combine externally.
	
  dput(list(coef=coef.triad1),"table1-ergm-ivnetwork1.txt")
  dput(list(coef=coef.triad2),"table1-ergm-ivnetwork2.txt")
  dput(list(coef=coef.triad3),"table1-ergm-ivnetwork3.txt")
  dput(list(coef=coef.triad4),"table1-ergm-ivnetwork4.txt")
  dput(list(coef=coef.triad5),"table1-ergm-ivnetwork5.txt")

# Now run the models with the instrumented score.

coef.triad1 <- NULL
coef.triad2 <- NULL
coef.triad3 <- NULL
coef.triad4 <- NULL
coef.triad5 <- NULL

m <- 1000
set.seed(10)
uyr <- unique(maoz.dat$year)
for(i in 1:1000){
  
	yrsi <- sample(uyr,length(uyr),rep=T)
	yrsil <- as.list(yrsi)
	inds <- unlist(lapply(yrsil,whichab,y=maoz.dat$year))
	datai <- maoz.dat[inds,]

	mod.triad1 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltrade_seq_hat1+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)
	mod.triad2 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltrade_seq_hat2+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)
	mod.triad3 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltrade_seq_hat3+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)
	mod.triad4 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltrade_seq_hat4+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)
	mod.triad5 <- glm(dichmid20~Lminreg302+Lcaprat301+Ldistance+Latopsecorr+Ligosecorr+Ltrade_seq_hat5+t1+t2+t3+peaceyrs+tchange+tschange, data=datai,family=binomial)

	coef.triad1 <- rbind(coef.triad1,rbind(mod.triad1$coefficients))
	coef.triad2 <- rbind(coef.triad2,rbind(mod.triad2$coefficients))
	coef.triad3 <- rbind(coef.triad3,rbind(mod.triad3$coefficients))
	coef.triad4 <- rbind(coef.triad4,rbind(mod.triad4$coefficients))
	coef.triad5 <- rbind(coef.triad5,rbind(mod.triad5$coefficients))

	}
	
# Output the results to access and combine externally.
	
dput(list(coef=coef.triad1),"table1-ergm-ivscore1.txt")
dput(list(coef=coef.triad2),"table1-ergm-ivscore2.txt")
dput(list(coef=coef.triad3),"table1-ergm-ivscore3.txt")
dput(list(coef=coef.triad4),"table1-ergm-ivscore4.txt")
dput(list(coef=coef.triad5),"table1-ergm-ivscore5.txt")
